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Abstract 

The problem of pattern selection arises when the evolution equations have 
many solutions, whereas observed patterns constitute a much more restricted 
set. An approach is advanced for treating the problem of pattern selection by 
defining the probability distribution of patterns. Then the most probable pattern 
naturally corresponds to the largest probability weight. This approach provides 
the ordering principle for the multiplicity of solutions explaining why some of 
them are more preferable than other. The approach is applied to solving the 
problem of turbulent photon filamentation in resonant media. 
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1 Introduction 



Nonequilibrium phenomena in real systems are usually described by sets of nonlinear 
differential or integro-differential equations in partial derivatives. The solutions to 
such equations are in many cases nonuniform in space exhibiting the formation of 
various spatial structures. It often happens that a given set of equations possesses 
several solutions corresponding to different spatial patterns. And, moreover, it is often 
impossible to distinguish between these solutions being based on stability analysis since 
all these solutions can be stable. If an ensemble of solutions can be parametrized by a 
multiparameter (5 from a manifold £>, so that all solutions corresponding to any f3 G B 
are stable, then the manifold B is called the stability balloon [1]. 

In the case of such a nonuniqueness of solutions one may wander which of the 
multiplicity of allowed states would actually be found for a given experimental protocol. 
This is how the problem of pattern selection arises, when the considered equations have 
many solutions for given external conditions, whereas observed patterns may constitute 
a more restricted set. And there no doubts are numerous systems and phenomena in 
real life with such a complicated behaviour [1]. 

Not only the solutions to equations can be nonunique but real phenomena may 
also display a variety of patterns for the same given conditions. This means that the 
multiplicity of solutions is not just an artifact but a feature expressing the intrinsic 
complexity of natural phenomena. If real patterns show multiplicity, is there any or- 
dering between them, such that one solution is preferred over the other? In equilibrium 
thermodynamics this is a familiar concept, with the free energy providing the ordering 
principle. But in nonequilibrium systems, there is no such a general organizing princi- 
ple to apply [1]. For some particular cases several heuristic recipes have been suggested. 
For example, one state may be preferred over the other if it has a larger basin of at- 
traction for typical initial conditions. Or one may think that the fastest growing mode 
dominates the evolution. Or selection is made assuming that the pattern with the 
fastest spatial decay rate is preferable. The limits of such heuristic arguments are dis- 
cussed in detail in Cross and Hehenberg [1] where one can find extensive citations to 
literature on the problem. 

In the present paper a general approach is developed for treating the problem of 
pattern selection. The approach, which is formulated in Sec. 2, is based on defining the 
probability distribution of patterns. In Sections 4 to 5, this approach is applied to the 
problem of turbulent photon filamentation. An introduction to the latter problem is 
given in Section 3. The obtained results are in very good agreement with experiment. 
The final Section 7 contains conclusions. 

2 Probability Distribution of Patterns 

Let us consider a system of evolution equations, which displays the multiplicity of 
solutions describing different spatial structures. Assume that these solutions can be 
parametrized by a multiparameter f3. All admissible values of (3 form a manifold 
B = {[3}. In many cases, differential equations in partial derivatives can be reduced to 
a d- dimensional system of ordinary differential equations, with a dimension d that may 
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equal infinity [2]. For the simplicity of notation, we shall keep in mind this possibility 
of working with a d-dimensional dynamical system. A generalization of this case will 
be given at the end of the present Section. 

The state of a rf-dimensional dynamical system is the set 



y(t) = { yi (t)= yi (l3,t)\i = 1,2,..., d} 



(1) 



of solutions to the system of differential equations, which can always be presented in 
the normal form 

j t y(t) = v(y,t), (2) 

where v = {vi\ i — 1, 2, . . . , d} is a velocity field. By assumption, the multiplicity 
of solutions is parametrized by a multiparameter (3, but in intermediate calculations 
we shall not, for the compactness of notation, label functions with f3, restoring this 
dependence in final formulas. 

Our aim is to classify the states (1) labelled by (3 by defining a probability measure 
on the manifold B. To introduce the probability distribution p{(3, t) of patterns labelled 
by a multiparameter (3, we may resort to the ideas of statistical mechanics [3], where 
a probability p can be connected with entropy S by the relation p ~ e~ s . In the 
nonequilibrium case, the entropy is a function of time, S(t). Since it is not the entropy 
itself but rather its change that is measurable, it is convenient to count entropy from 
its initial value S(0), thus, considering the entropy variation 



AS(t) = S(t) - S(0) , 



(3) 



which is a kind of relative entropy [4]. The entropy may be defined as the logarithm 
of an elementary phase volume [3], which for the nonequilibrium case can be expressed 
as 

S(f)=ln|*T(f)|, (4) 



with the elementary phase volume 

5r(t) = H5 yi (t). 

i 

Hence the entropy variation (3) is 

AS(t) = In 



(5) 



ST{t) 



5T(0) 



(6) 



AS 



The probability distribution p ~ e 
form 

e -A5(/3,t) 



, being normalized with respect to (3, takes the 



p((3,t) 



Z(t) 



Z(t) 



dp , 



(7) 



where the integration over (3 runs through the manifold B. If the latter manifold is not 
continuous, the integration is to be replaced by summation. 
The elementary phase volume (5) can be presented as 



5r(t) = nE M y -(%(°). 



(8) 
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where My(t) are the elements of the multiplier matrix [5] defined through the varia- 
tional derivatives 

M^) = M|, A-/ y (0) - S u . (9) 
Then the entropy variation (6) writes 

AS-(t)=I>l M «(*)l- (10) 

i 

The multiplier matrix M(t) = [Mij(t)\ satisfies the equation 

j t M(t) = J(y,t)M(t), (11) 

which follows from the variation of the evolution equation (2) and where J(y, t) = 
[Jij(y,t)] is the Jacobian matrix with the elements 

j , A Svj(y,t) 

Jy(y '* ) = "W (12) 

The initial condition for Eq. (11), according to the definition (9), is M(0) = [Sij]. 
Combining Eqs. (7) and (10), we get the probability distribution of patterns, 

*^nRM' (13) 

with the normalization factor 

z(t)= /n \M^M\ dfl - 

This is the probability distribution of solutions classified with a multiparameter (3. 
Since each solution, by definition, represents a particular pattern, expression (13) is 
the probability distribution of patterns. This expression naturally connects the notion 
of probability and the notion of stability. Really, the multipliers are smaller by modulus 
for more stable solutions and, consequently, for more probable patterns. 

To calculate the probability distribution of patterns (13), we need to know the 
multipliers (9) which are defined by the equation (11). Using the latter equation, 
the pattern distribution (13) can be transformed as follows. Introduce the matrix 
L(t) = [Lij(t)\ with the elements 

Ly(t)=]n|JUy(t)|. ( 14 ) 
Then the entropy variation (10) writes 

AS(t) = Tr L(t) . (15) 

The trace of a matrix does not depend on the matrix representation. Hence, we may 
perform intermediate calculations using one particular representation, returning at 
the end to the form independent of a representation. To this end, let us consider a 
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representation when the multiplier matrix is diagonal. Because of Eq. (11) with the 
initial condition My(0) = 5ij, the multiplier matrix is diagonal if an only if the Jacobian 
matrix is diagonal too. In this case, the evolution equation (11) yields 

M u {t) = exp y* Ju(y(t'),t') dt'} . 

From Eq. (14) we have 

Lii(t) = f ReJ u (y(t'),t') dt' . 
Jo 

Introducing the notation 

tf(t) = E ReMy,t), (16) 

i 

we get 



Tr Lit) = ( Kit') dt' 
Jo 



Without the loss of generality, we may assume that the state (1) consists of real func- 
tions since any complex function can always be treated as a pair of real functions. 
Hence the velocity field in Eq. (2) can also be considered as real. Then the eigenvalues 
of the Jacobian matrix (12) are either real or, if complex, come in complex conjugate 
pairs. Therefore 

E Re My, t) = E My, t) = Tr J( y , t) . 

i i 

Thus, Eq. (16) can be written as 

K(t)=TrJ(y,t). (17) 
And for the entropy variation (15), we find 

ASit) = C Kit!) dt! . (18) 
Jo 

Expression (17) in dynamical theory is called the contraction rate [6]. All consider- 
ation given above can be straightforwardly generalized to the case when the state (1) 
consists of functions yi(x, t), where x is the set of space variables. This results not more 
than in a slight complication of notations. Then everywhere the variable x appears 
together with the index i as an additional continuous index, and the sums over i are to 
be complimented by the integrals over x. The multiplier matrix (9) becomes a matrix 
with respect to indices i, j as well as with respect to x, x', 



with the initial condition 

Mij(x, x', 0) = Sij5(x — x') . 
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Similarly, the Jacobian matrix (12) becomes a matrix with the elements 

Jij{x,x ,y,i) — . . . 

oyj{x',t) 

Employing the matrix notation with respect to continuous variables [5] , we may repeat 
the same steps as above in deriving the probability distribution (13). As is mentioned, 
the sum over x is to be understood as the corresponding integral. And the product 
over a continuous variable can be defined [7] as 

jj [ f(x) = exp / In f{x) dx . 
As a result, the pattern distribution (13) reads 

p (^'*) = exp {~51 / ^\Mu(x,x,(3,t)\ cfoj . 

The contraction rate (17) becomes 

K(t) = J2 I Mx,x,y,t)dx, (19) 
i J 

where y = {yi(x, t)}. In this way, the contraction rate has always the form of Tr J(y, t), 
as in Eq. (17), where the trace has to be defined according to the representation of the 
Jacobian matrix. 

Restoring in the contraction rate the dependence on the parameter f3 labelling 
different patterns and using the entropy variation (18), we finally obtain the probability 
distribution of patterns 

p(M = ^ ex p{-^ K (PJ) dt '} > (2°) 

in which the normalization factor is 

Z(t) = J exp j-^* K((3,t') dt'} d/3 . 

Thus, each solution labelled by (5 is equipped with the probability weight (20). 
Consequently, that pattern is preferred over the other which has a higher probability 
weight. This is equivalent, because of the form (20), to saying that one pattern is 
preferable over others if its local contraction 

A((3,t)= i f K((3,t')dt' 
t Jo 

is minimal with respect to (3. The latter provides the ordering principle for classifying 
solutions and related patterns. The local contraction plays for nonequilibrium dynam- 
ical systems the role analogous to that of the free energy for equilibrium statistical 
systems. 
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3 Turbulent Photon Filamentation 

To illustrate the probabilistic approach to pattern selection, developed in the previous 
section, let us consider spatial structures appearing in resonant samples when increas- 
ing the Fresnel number. In active media interacting with electromagnetic field, there 
appears a variety of spatiotemporal structures. For example, electric field in laser 
cavities can exhibit spatial states, such as solitons and vortices, which bear a close 
analogy with similar structures in liquids [8,9]. In general, there is a direct correspon- 
dence between the Maxwell-Bloch equations for slowly varying field amplitudes and 
hydrodynamic equations for compressible viscous liquid [1,10]. The Fresnel number for 
optical systems plays the same role as the Reynolds number for fluids. In the same 
way as when increasing the Reynolds number, the fluid becomes turbulent, the field 
dynamics exhibits chaotic behaviour when increasing the Fresnel number. Similarly to 
hydrodynamics, optical systems display spatial multi-stability with several coexisting 
distinct stable states for the same values of parameters [1,11]. Spatiotemporal chaos 
in optical systems is characterized by the same fast decay of correlation functions as 
it occurs for turbulent fluids. This is why one calls such chaotic optical phenomena 
optical turbulence [1,10,11]. 

In active optical systems, having the standard cylindrical shape, the transition from 
regular behaviour to spatiotemporal chaos occurs with the increasing Fresnel number 




where R is the cylinder radius, L is the characteristic length, and A is the radia- 
tion wavelength. Physical processes accompanying the route from a regular regime to 
chaotic one are similar in different active optical media. 

At small Fresnel numbers F < 1, there exists the sole transverse central mode 
practically uniformly filling the medium. When the Fresnel number is around F ~ 1, 
the cavity can house several transverse modes seen as a regular arrangement of bright 
spots in the transverse cross-section. These regular spatial structures emerge from 
an initially homogeneous state with a break of space-translational symmetry. They 
are regular in space forming ordered geometric arrays, such a polygons, and they are 
regular in time being either stationary or periodically oscillating. These transverse 
structures are imposed by the cavity geometry and correspond to the empty-cavity 
Gauss-Laguerre modes. Such regular structures are well understood theoretically, their 
description being based on field expansions over the modal Gauss-Laguerre functions 
prescribed by the cylindrical geometry, and the theory being in reasonable agreement 
with experiments for lasers, e.g. for CO2 and Na2 lasers [11-17], and for active nonlinear 
media, as the photorefractive Bi 12 Si02o crystal pumped by a laser [18-20]. Similar 
structures also arise in many passive nonlinear media, such as a Kerr medium [20]. For 
Fresnel numbers up to F ~ 5, the number of bright spots is proportional to F 2 . In the 
longitudinal cross-section this corresponds to the existence of bright filaments whose 
number follows the F 2 law as F increases. 

Around F m 10 there occurs a principal change of properties, from the existence 
of regular structures to a turbulent-type state [21-23], with the intermittent behaviour 
in the region 5 < F < 15. This happens in lasers [21-23] as well as in photorefractive 
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crystals [18-20]. For Fresnel numbers F > 15, the arising spatial structures are very 
different from those associated with the empty-cavity modes. The modal expansion is 
no longer relevant and the boundary conditions have no importance. The medium looks 
as consisting of a large number of parallel independently flashing filaments, whose num- 
ber is proportional to F, contrary to the case of small Fresnel numbers with the number 
of filaments proportional to F 2 . The filaments are chaotically distributed in space, are 
not correlated with each other, and are aperiodically flashing in time. Such a spatio- 
temporal chaotic behaviour is characteristic of hydrodynamic turbulence, this is why 
the similar phenomenon in optics is commonly called the optical turbulence [1,10,11,23]. 
Contrary to the regime of small F, when the regular spatial structures are prescribed 
by the geometry and boundary conditions imposing their symmetry constraints, the 
turbulent optical filamentation is strictly self-organized, with its organization emerging 
from intrinsic properties of the medium [20,24]. This kind of optical turbulence has 
been observed in both photorefractive crystals [18-20] and lasers [21-23,25-30]. Espe- 
cially accurate and thorough experimental studies for CO2 and Dye lasers have been 
accomplished in Refs. [25-30]. The independence of turbulent optical filamentation of 
boundary conditions and, hence, its purely self-organized nature are confirmed by the 
observation of this effect in the resonatorless discharge-tube superluminescent samples, 
such as lasers on Ne, Tl, Pb, N 2 , and N2" vapors [31-35]. Since the optical turbulence 
is characterized by the formation of bright filaments with a high density of photons, 
this phenomenon may be named the turbulent photon filamentation. 

Let us summarize the main features characterizing the regular and turbulent regimes 
in active optical media. At low Fresnel numbers F < 10, the regular regime occurs 
whose typical characteristics are: (i) Bright filaments are regularly arranged in space 
forming an ordered structure seen in the transverse cross-section as a polygon made 
of bright spots, (ii) The spatially ordered structure is either stationary or periodically 
oscillating in time, (iii) The number of bright filaments is proportional to F 2 . 

At high Fresnel number F > 10, the turbulent regime develops whose typical fea- 
tures are: (i) The bright photon filaments are chaotically distributed in space, (ii) The 
filaments aperiodically flash in time, (iii) The number of these uncorrelated filaments 
is proportional to F. 

While for the low-Fresnel-number regime, with regular optical structures, a good 
theoretical understanding was achieved [11-20], for the high-Fresnel-number regime of 
the turbulent photon filamentation, there has been no persuasive theory developed. 
This problem was addresses in Refs. [36-39], where the consideration was based on a 
simple model, only the stationary case was analysed, and the minimum-energy argu- 
ments were employed. This model consideration showed that photon filaments can re- 
ally be formed in resonant media due to an effective interaction between atoms through 
photon exchange. The estimates for the filament radius turned out to be in good agree- 
ment with experiment. However, the empirical approach of Refs. [36-39] cannot be 
accepted as completely satisfactory. This is mainly because of the two principal points. 
One is that a stationary case was addressed, while turbulence is rather a nonstationary 
phenomenon and must be treated by time-dependent evolution equations. The sec- 
ond point is that the minimal-energy condition was invoked, while there is no such a 
principle for nonequilibrium systems. 

Here we demonstrate that the problem of turbulent photon filamentation can be 
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successfully treated by the probabilistic approach to pattern selection formulated in 
Section 2. Our consideration is based on realistic evolution equations for resonant 
media. Such equations are necessary for calculating the contraction rate K defining 
the probability distribution of patterns (20). 



4 Equations for Resonant Medium 

The system of iV resonant atoms interacting with electromagnetic field is described by 
the Hamiltonian 

H = H a + H f + H af , (21) 
consisting of the following terms: The Hamiltonian of two-level resonant atoms 

1 N 

H a = ^Y. wb(l + 0, (22) 

with the transition frequency u , where of is a Pauli operator. The radiation-field 
Hamiltonian 

Hf= ^ / {E 2 + H") dr, (23) 

— » — * — * — * — * 

with electric field E and magnetic field H = V x A, where the vector potential A is 
assumed to satisfy the Coulomb gauge calibration V- A = 0. The atom-field interaction 
Hamiltonian 

H af = - E (-jj-Ai + di-Eoi) (24) 

— * — * 

has the standard dipole form, in which A t = A(fi,t); the transition-current and tran- 
sition dipole operators are 

ji = iooq (daf — d*a~^j , di = daf + d*a~ , 

where d is the transition dipole and of are the rising or lowering operators, respectively; 
and E 0i is a seed field [40]. 

To consider nonuniform systems with arising spatial structures, it is convenient to 
invoke the space-time representation of evolution equations, whose general description 
is given in the book [40] and a detailed analysis in Refs. [41,42]. The equations are 
written for the average quantities 

u(f,t) = < a~(f,t) > , s(f,t) = < cr z (f,t) > , (25) 

using the standard semiclassical and Born approximations. For the compactness of 
presentation, let us introduce the notation 

f(f,t) = f (f,t) + f rad (f,t) (26) 

for an effective field acting on an atom. This field consists of the term 

fo(f,t)= -id-E (r,t), (27) 
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due to the cavity seed field, and of the term 



frad(r,t)= - \%1P J [p{f-f')u{f',t)-e d 2 ^{f-f')u*{f',t)\ dr> (28) 

corresponding to the interaction of atoms through the common radiation field, where 
p is the density of atoms, d = doed, do = \d\, and 

expgfep \f\) k -^o 4 3 2 

The seed field 

E (f, t) = E x e i{kz ~ wt) + E* e"*^-^) (29) 

selects a longitudinal mode propagating along the axis z which is the axis of a cylindrical 
sample. The resulting equations are 

du . 

— = -{lUjQ + 7 2 )m + sf , 

ds 

J= -2(«V + /*«)-7i(a-0, (30) 

^= -2 72 H 2 + S K/ + r«), 

where 71 and 72 are the level and line widths, respectively, and ( is the pumping 
parameter defined by the intensity of nonresonant lamp pumping. 

The sample has cylindrical shape typical of lasers. The radiation wavelength A, 
cylinder radius R, and its length L are related by the inequalities 

I « 1 , I « 1 ■ (31) 
We consider the quasiresonance situation when 

^ < 1 , A = to - uj . (32) 

As usual, the relaxation parameters are assumed to be small, so that 

— < 1 , — < 1 . (33) 

The solutions to Eqs. (30) are, in general, nonuniform. This is because atoms 
interact with each other thorough the radiation field (28) containing the function <f(f) 
fastly oscillating in space and diminishing with increasing |f|. Hence the atoms, situ- 
ated far from each other, do not interact and radiate independently, although in the 
longitudinal direction they may be correlated due to the cylindrical geometry of the 
sample and the propagation of radiated field along the axis z. This suggests to look for 
a solution of Eqs. (30) in the form of a bunch of filaments aligned along the cylinder 
axis. Let us imagine that the radiation field inside the sample is stratified into Nf pho- 
ton filaments stretched along the axis z. These filaments can have different radii r n , 
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with n — 1,2, ... , Nf. The radiation is mainly concentrated inside the bright filaments, 
fading away outside them, so that at the radial distance b n the corresponding solutions 
are an order of magnitude smaller than at the axis of a filament. The relation between 
the characteristic lengths b n and r n can be written down if the profile of solutions inside 
a filament is known. If the filament profile can be approximated by the normal law 
exp(— r 2 /2r^), with the filament radius r n being the standard deviation, then 



b n = V2 In 10 r n . (34) 

The filaments do not interact with each other, because of which their locations in the 
radial cross-section are random. This means that the filament axes are located at 
random points {x n , y n }, where n — 1,2, ... , Nf. 

The above discussion explains why it is reasonable to try to find the solutions of 
Eqs. (30) in the form of expansions 



N f N f 

u(r, t) = J2 «n(r, t)Q n (x, y) e lkz , s(r, t) = £ s n (r, t)Q n (x, y) (35) 

n=l n=l 

over filaments, where 

G n (x, y) = G (b 2 n -{x- x n f + (y- y n f) 

is a unit-step function. Substituting the presentation (35) in Eqs. (30), we obtain a 
system of equations for u n , s n , and |w n | 2 . These equations can be simplified if we are 
not interested in the detailed internal structure of each filament but rather wish to find 
out their characteristic sizes. Then we may employ the mean-field approximation for 
the averages 

u(t)= ^- I u n (f,t)df, s(t)= ^- f s n (f,t)df, 

\u(t)\ 2 = ±r ( \u n {r,t)\ 2 df, (36) 



V n JVr 

where the integration is over the volume V n = irb n L and the index n in the left-hand 
sides, for short, is dropped. 

To present the equations for the functions (36) in a compact form, we introduce 
two effective coupling parameters: 

9= -Prr I "~ L " UI ' , ,V :r ~ n dr&' (37) 



472 V n Jv, 

and 

!J 



sin[/c 




- k(z - z')] 




k Q \f - 




cos[/cq 




- k(z - z')\ 




k \f- 


■r'\ 



These parameters enter the definitions of the collective frequency and collective width, 
respectively, 

ft = wo + $V , r = 72 (1-^). (39) 
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Then from Eqs. (30), employing the introduced notation, for the functions (36) we 
obtain the equations 

^ = -m + T)u - isd- E x e' wt , 
dt 

" S -Ag l2 \u\ 2 - ^(s - C) - 4 Im (u*tf • .Ei e**) , (40) 



dt 

d\u 



-2Y\u\ 2 + 2s Im [u*d- E x e~ iult ) , 



dt 

describing the time evolution of a filament. Such a reduction of initial equations for 
a bunch of filaments to the set (40) of equations for each of the filaments has become 
possible due to the absence of interactions between filaments. Different filaments can 
have different radii whose distribution is yet to be defined. 

It is worth emphasizing that the presentation of solutions as the expansions (35) 
over filaments is rather general and includes as well the case when there are no separate 
filaments at all. Really, if it turns out that the most probable filament radius is close 
to the radius of the sample, R, this would actually mean that there are no several 
filaments, but the whole volume of the sample is filled by one filament. 

5 Characteristics of Arising Filaments 

To analyze Eqs. (40), we employ the scale separation approach [43,44] which is a 
generalization of the averaging technique [45] to the case of nonequilibrium statistical 
systems. To this end, we notice first of all that the functional variables in Eqs. (40) can 
be classified onto fast and slow in time. Due to the existence of the small parameters 
(33), the variable u is fast as compared to the slow variables s and \u\ 2 . The slow 
variables play the role of quasi-invariants for the fast function u. With s being a 
quasi- invariant, the solution for the fast variable u, described by the first of Eqs. (40), 
writes 

— * — * 

u(t) = u e -^+ r )* + \e- wt - e -^+ r )'j , (41) 

uj — \ l + u L J 

where u = u(0). This solution is to be substituted in the second and third of Eqs. 
(40) whose right-hand sides are to be averaged over time explicitly entering the fastly 
oscillating functions. In this way, we meet the quantity 

Im 1 r T -> -» 

a = — lim - / u*(t)d- E x e' luJt dt (42) 



sT T ^°° T Jo 

describing the influence of the seed field on atoms. Taking into account Eq. (41), we 
have 

\d-E!\ 2 

a= (J-ar + T* - (43) 

Recall that the seed field is, by definition, a weak cavity field whose role is to fix the 
axis of the cylindrical sample. This field does not pump the energy into the system 
which is standardly done by nonresonant lamps and which is described by the pumping 
parameter ( in the evolution equations (30) or (40). The seed field selects a longitudinal 
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mode imposing no restrictions on transverse modes. The selection of the latter has to 
be done by the internal properties of the evolution equations. The weakness of the 
seed field is to be understood in the sense of the smallness of the parameter (43), 

a < 1 . (44) 

In analysing further the evolution equations (40), it is convenient to make the 
transformation of the functional variable \u\ 2 to 

w = \u\ 2 — as 2 . (45) 

Finally, following the scale separation approach [43,44], from the last two of Eqs. (40) 
we obtain the equations 

ds , dw n . , 

_ = -Ag^ 2 w -7i(s - C) , — = -2-f 2 {l - gs)w (46) 

for the slow variables. These equations can be written in the form (2), 

ds dw 

Tt =Vl > Tt =v ^ 

with the velocity field 

vi = -Ag-f 2 w - 7i(s - C) , v 2 = -272(1 - gs)w . 

For the Jacobian matrix (12), we now have 

dv x dvi 
Ju = -^- = -7i , Jn = ^— = -4fi-72 , 
OS ow 

dv 2 dv 2 
■hi = -7^7 = -272W , J 22 = = -272(1 - gs) . 

From here, we find the contraction rate (17), 

K = Tr J = - 7l - 2 72 (1 - gs) . (47) 

The value of Eq. (47) depends on the characteristics of filaments through the coupling 
parameter (37) containing in its definition the volume V n = ^h 2 n L, which is the volume 
of a cylinder enveloping a filament. The radius of the enveloping cylinder, b n , is related 
to the filament radius, r n , by relation (34). For what follows, it is convenient to 
introduce the dimensionless quantity 

Thus, the contraction rate (47) is a function K = K((3,t) of the variable (48) and also 
a function of time entering through s = s(t). For different filaments, the variable (3 
can take different values from the interval 

< < F , (49) 
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where F is the Fresnel number. The distribution of the variable (3 in the interval (49) 
is given by the probability distribution (20). The maximum of the latter defines the 
most probable f3 and, respectively, the most probable filament radius. 

Let us consider the very beginning of the process, when t — > +0, in order to 
understand what are the most probable filaments to be formed. At this initial stage, 
the probability distribution (20) can be written as 

P(M- exp{-K((3,0)t} . 

Hence the maximum of this distribution corresponds to the minimal contraction rate. 
The latter, according to Eq. (47), depends on (3 through g = g((3). We consider the 
standard laser-type setup, when the atoms at the initial time are not inverted, i.e. 
So < 0, and the medium is pumped with a long nonresonant pulse whose action is 
described by the pumping parameter ( > in the evolution equations (30), (40), and 
(46). In this case, the minimum of the contraction rate K(/3,0) is equivalent to the 
maximum of the coupling g(f3), which is defined by the equations 

dq d 2 q , . 

i= ' «! <0 - (50) 

To solve these equations, we need to analyse the integral (37) giving g = g(/3). 
The coupling g(fl) defined by the integral (37) can be presented in the form 



37T 7 pL 
S(J3) = TnB 



tt/3 — / Si(x) dx 
Jo 



(51) 



The procedure of reducing Eq. (37) to the form (51) is described in the Appendix. 
The extrema of g((3) are given by the equation 

Si(2/3) = 7i . (52) 

From several solutions of Eq. (52) we have to choose the absolute maximum of g(f3). 
The latter occurs at f3 — 0.96, which, according to Eq. (48), gives b n = 0.55\^XL. 
Using the relation (34), we find the most probable radius of a filament 



r f = 0.26VAL . (53) 

In the typical experiments [21-23,25-35] observing the turbulent photon filamentation, 
the excitation was achieved by means of the quasistationary nonresonant pumping. 
Such a pumping can be treated as a quasistationary process if its duration is much 
longer than the characteristic time 2h/uj of fast oscillations, which is always the case. 
The quasistationary pumping can be characterized by an effective pumping parameter 
( in the evolution equations. To describe the finite duration time of the pumping 
procedure, one can consider ( = ((t) as a slow function of time, slow in the sense of 
the inequality 

2tt d( 

— -r < 1 . 

uj dt 
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In this way, for the standard laser setup, the filament radius (53) defines the most 
probable pattern for the arising bunch of photon filaments. 

The most probable number of filaments can be evaluated from the normalization 
condition 

r s(r,t)dr = ((t), (54) 



V 

where the integration runs over the whole volume of the sample, V = nR 2 L. To 
approximately calculate the integral (54), let us consider the moment of time when the 
filaments have already been formed and the population difference inside each filament 
of the radius 77 has reached a value close to the pumping parameter (. Then we may 
write 

' s(f, t) dr~ N f V f C +{V- N f V f )s out , 



where Vf = nr 2 L and s out is a population difference outside the filaments. In this case, 
independently of the value s out , condition (54) yields 




N f =[-\ . (55) 



For the filaments of the most probable radius (53), equation (55) gives 

N f = 4.71F . 

The number of filaments is proportional to the Fresnel number, which is in complete 
agreement with all experiments on the turbulent photon filamentation. The coefficient 
4.71 is also in good agreement with experiments [18-20]. 



6 Dynamics of Filament Flashing 

The time evolution of the filaments is described by Eqs. (46). At the initial stage of 
the process, when 71 £ <C 1, we may omit the term containing 71. Then the resulting 
system of nonlinear equations can be solved exactly yielding the solutions 



— tanh ( -~\ + - , 



w 



7 ° sech 2/t - t0 



#72 V t J g ' 4# 2 7f V - (l 

in which 7 is the radiation width, r is the radiation time, so that 

7o 2 = T 2 + Ag 2 ^ (\u \ 2 - a sl) , 

where uo = u(0), so = s(0), ao = a(0), and 

ro = 72(1 -gs ) , 7oTo = i; 

the delay time to being given by the expression 



t n = — m 



70 - r 
70 + r 



(56) 



(57) 



(58) 



(59) 
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Invoking the transformation (45), we get 



| u |2 = Jfi S ech 2 ( t —^ ) + as 2 . (CO) 



The found solutions depend on the coupling g both directly as well as through the 
characteristic parameters 70, r , and to- This means that the properties of filaments 
with different radii are essentially different, since g = g((3). 

At the later stage, the term with 71 cannot be dropped. Then Eqs. (46) cannot 
be solved exactly. But we may look at the long-time behaviour, when the system 
approaches a stationary regime. The latter can be achieved if the duration of the 
pumping pulse is longer than 7^ as well as 7^ , which is usually the case when pumping 
a resonant medium. 

Before the stationary regime is reached, the oscillation of solutions is described by 
the eigenvalues of the Jacobian matrix associated with Eqs. (46). These eigenvalues 
are 

A ± = - \ [71 + 72(1 - gs)\ T {[7i - 272(1 - gs)] 2 - S2 1 2 2 g 2 w} 1/2 . (61) 

Because of the dependence of s(t) and w(t) on time, Eq. (61) shows that A ± (t) is also 
time dependent, which means that the oscillation of solutions cannot be defined by 
fixed frequencies but occurs in an aperiodic way. 

Al long times, such that 71 £ ^> 1 and 7 2 t 3> 1, the system approaches a stationary 
regime. There are two types of stationary solutions to Eqs. (46), one of them is given 
by the fixed point 

sI = C, < = 0; (62) 



and another, by the fixed point 



" 2 = ~i^r ' (63) 



The stability of the stationary solutions can be determined from the Lyapunov analysis. 
To this end, we have to evaluate the Jacobian eigenvalues (61) at the corresponding 
fixed points (62) and (63), which results in the characteristic exponents 



A+ = - 7 i, K = -272(1-00 (64) 



and, respectively, 




A 2 ± 

7i 




1 + 8^(1-00 >■ (65) 



The real parts of Eqs. (64) and (65) define the Lyapunov exponents whose signs charac- 
terize the stability of the related fixed points. As far as the values of the characteristic 
exponents (64) and (65) depend on the coupling g = g([3) which, according to Eq. 
(51), essentially depends on the parameter f3, that is, on the radius of a filament, the 
stability properties for the filaments of different radii can be drastically different. For 
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the parameter (3 in the interval (49), the coupling (51) varies in a wide diapason. Thus, 
for small (3, we have 

I672 

hence g((3) — > as (3 — > 0. For large /5, Eq. (51) gives 

In the case of the wavelength A>a, much larger than the mean interatomic distance, 
one has pA 3 1. Since usually 7 ~ 72 and L ^> A, the coupling may reach quite 
large values g((3) 3> 1. Such a wide variation of the coupling g = g(f3) results in 
rather different characteristic exponents (64) and (65) for different filaments and, as a 
consequence, in essentially different stability properties of the latter. 

All filaments can be separated into three types. One group consists of those fila- 
ments for which 

g(P)C < 1 . (66) 

Then the fixed point (62) is a stable node, while that (63) is a saddle point. The 
filaments with the radii satisfying condition (66) are characterized by the stationary 
solutions (62), which shows that these filaments after the time 7\ = 7f 1 practically 
stop radiating coherently. 

Increasing g(/3), one reaches the equality g{P)C = 1, when both fixed points (62) 
and (63) merge together becoming neutral. In the interval 

1 < g{m < 1 + , (67) 

072 

the fixed point (62) is a saddle point and that (63) is a stable node. Hence, the filaments 
with the radii satisfying condition (67) are described by the stationary solutions (63). 
Such filaments continue radiating coherently in the stationary regime, although the 
level of coherence, characterized by the value of u^, is not high. 

The third group of filaments corresponds to large coupling, such that 

9(P)( >!+«£-■ (68) 
072 

Then the fixed point (62) remains a saddle point while that (63) becomes a stable focus. 
The latter means that the approach of the related solutions to the stationary point (63) 
is not monotonic but through a series of pulses. The filaments with the radii satisfying 
condition (68) radiate by bright flashes interrupted by the intervals of darkness. This 
intermittent behaviour, for finite times, is not periodic, as is seen from the form of Eq. 
(61), but as t — > 00, there appears an asymptotic period. This can be noticed from the 
characteristic exponent (65) which, for the case (68), can be written as 

Af = - y =F iuco , (69) 



where the asymptotic frequency is 




1/2 



T 1 8 ^ b(/3)C - 1] - 1 j • (70) 
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The value of the latter, because of the dependence on <?(/?), is essentially different for 
different filaments. 

7 Conclusion 

An approach is developed for treating the problem of pattern selection. The approach 
is based on defining the probability distribution of patterns. This gives the ordering 
principle for the multiplicity of admissible solutions to nonlinear differential equations 
and, respectively, to the corresponding space structures. The maximum of the proba- 
bility distribution naturally defines the most probable pattern. 

This probabilistic approach is applied to the problem of turbulent photon filamen- 
tation in laser media. The most probable filament radius is found, and the number of 
filaments is evaluated. The found values are in very good agreement with experiment. 
Thus, the filament radii observed in experiments with dye and CO2 lasers [25-30] are 
77 ~ O.Of cm for dye lasers and 77 ps 0.08 cm for CO2 lasers, which is in perfect 
agreement with formula (53). The filament radii measured in experiments with Ne, Tl, 
Pb, N 2 , and N 2 vapor lasers [31-35] are 77 w 0.01 cm, which again coincides with the 
result of Eq. (53) for the corresponding wavelengths and laser lengths. The number of 
filaments, given by Eq. (55) is Nf ps 10 2 — 10 3 for different experiments [31-35]. The 
found dependence of the filament radius 77 ~ F" 1 / 2 and filament number Nf ~ F on 
the Fresnel number is the same as observed experimentally. 

The turbulent photon filamentation provides a good example of how a rather com- 
plicated phenomenon can be successfully described by the probabilistic approach to 
pattern selection. The formulation of this approach in Section 2 is general, which per- 
mits one to apply it for different phenomena where the problem of pattern selection 
arises. 
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Appendix. Effective Coupling 

The effective coupling g is given by the integral (37) over the enveloping cylinder 
of radius b n . Assuming that A <C b n , the integral (37) can be reduced to the form 



37T7p f b " , f L / 2 sm(k y/r 2 + z 2 - kz) , 

g = / r dr dz , 

272 J J-l/2 koVr 2 + z 2 

where r is the radial variable. Because of the quasiresonance condition (32), we have 
k w k. Then, with the change of the variable x = k(y/r 2 + z 2 — z), we get 

37T7p f b ™ r kL sin a; 
g = — / r ar dx . 

272K J Jkr 2 /L X 

Here, the integration limit kL, due to the inequality A <C L, can be replaced by 00. As 
a result, we have 



37T7P f 6 " 

9 = 



2 l2 k 
where 



/ 

Jo 



7T n . /Ar 2N 
2" Sl 



N f x sint , 7r p sint , 
Si(x) = / dt = - + / di 

y ' JO t 2 ioo t 

is the integral sine. Introducing notation (48), we come to g{(3) given in Eq. (51). In 
the same way, for the coupling (38) we find 



where 

r x cost , 



ca(x) = r 

JO 



t 

is the integral cosine. The found expressions for the effective couplings can also be 
transformed by means of the integrals 

J Si(x) dx = xSi(x) + cosx , J Ci(x) dx = xG\{x) — sinx . 

This yields 

9(P) = ^Pl^ W ~ W Si (2/9) + 1 - 008(2/3)] , 
472 Ar 

</(£) = [sin(2/3) - 2/5 Ci(2/5)] . 
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